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NOMENCLATURE 

A  =  amplitude  of  sinusoidal  oscillations 

a  =  transformation  parameter 

Cjl  =  inline  force  coefficient 

Cl  =  transverse  force  coefficient 

Cj  =  Fourier-averaged  drag  coefficient 


-m 


inertia  coefficient  =  1  +  k^ 


D  =  diameter 

Fjl  =  inline  force 

Fl  =  transverse  force 

K  =  Keulegan-Carpenter  Number  =  UmT/D 

K+  =  K  by  current  =  K(l  +  U0/Um)  =  K  +  U0T/D 

kV.  =  Fourier-averaged  added  mass  coefficient 

n  =  order  of  polynomial  (see  Eq  17) 

ps  =  pressure  on  cylinder 

Poo  =  pressure  at  outer  boundary 

ps  =  pressure  coefficient  =  (ps  -  p^,)  /  pUoo 

R  =  cylinder  radius 

Re  =  Reynolds  number  VU/v,  or  UmD/v 


IX 


Re+  =  Re  modified  by  current  =  Re(l  +  U0/Um)  =  Re  +  U0D/v 

r  =  radial  distance 

T  =  period  of  oscillations 

t  =  time 

U  =  time  dependent  velocity 

Um  =  maximum  velocity  in  pure  sinusoidal  flow 

Uo  =  collinear  steady  current  velocity 

V  =  constant  velocity  reached  at  the  end  of  the  acceleration  period 

Vr  =  current  ratio  =  U0/Um 

p  =D2/vT  =  Re/K 

A  a  =  disturbance  oscillation  applied  to  the  ambient  flow,  (in  degrees) 

Ac,  =  computational  grid  spacing 

Q  =  relative  displacement  of  fluid  during  the  acceleration  period 

=  (S/R)v  =  0.5  (dU/dt)  tJ/R  =  0.5  V  tv/(R) 

]i  =  dynamic  viscosity 

v  =  kinematic  viscosity 

p  =  density 

0  =  angular  position 

\\f  =  stream  function 

co  =  vorticity 


ACKNOWLEDGMENTS 

For  his  support  and  guidance  both  as  teacher  and  friend  I  express  my 
sincere  appreciation  to  Distinguished  Professor  T.  Sarpkaya. 

Many  others  have  given  me  help  along  the  way.  Dave  Marco  and  Tony 
Circelli  have  given  me  tremendous  support  using  the  ME  and  AERO 
computer  labs.  John  Ekaterinares  and  Dr.  Kwok  of  NASA  provided  the 
avenue  for  use  of  the  CRAY  supercomputers.  LT  K.  Charles  Vanover  was 
instrumental  in  my  introduction  to  the  existing  code.  Professor  Dalton  and 
Dr.  Xuegeng  Wang  supplied  the  initial  program. 

Finally  my  loving  appreciation  goes  to  my  wife  Jan  for  her  continual 
support  and  understanding. 


XI 


I.     INTRODUCTION 

Numerical  experimentation  in  fluid  dynamics,  through  the  use  of 
finite-difference,  finite-element,  and  discrete-vortex  methods,  has  attracted 
considerable  attention  during  the  past  two  decades  and  produced  laminar 
flows  difficult  to  measure  and  turbulent  flows  hard  to  verify  and  impossible 
to  generalize.  The  reasons  for  this  are  relatively  simple.  Numerical  solutions 
based  on  the  full  Navier-Stokes  equations  are  not  stable  at  high  Reynolds 
numbers  and  the  instability  is  non-linearly  related  to  the  particular  flow, 
input  parameters  and  the  discretization  conditions.  Also,  the  real  flow  at  the 
computed  Reynolds  numbers  may  be  turbulent,  at  least  in  some  regions  of 
the  flow,  and  the  numerical  experiment  does  not  imitate  the  physical 
experiment.  Furthermore,  the  observed  physical  and  numerical  instabilities 
do  not  necessarily  correspond  to  each  other.  Assuming  that  the  calculations 
for  a  given  flow  are  carried  out  at  sufficiently  small  Reynolds  numbers,  where 
the  flow  is  known  to  remain  stable  and  laminar,  one  quickly  discovers  that  it 
is  practically  impossible  to  measure,  to  any  credible  degree  of  accuracy,  most 
or  all  of  the  predicted  quantities  (except  the  Strouhal  number  and  the 
photographs  of  the  flow  patterns). 

Evidently,  one's  view  of  the  state  of  the  numerical  modelling  depends  to 
a  large  extent  on  one's  objectives.  For  example,  if  the  objective  is  to  obtain 
some  approximate  answers  and  flow  kinematics,  one  might  be  perfectly 
satisfied  with  the  existing  codes.  If  the  objective  is  to  match  the  measured  and 
calculated  results  (e.g.,  lift  and  drag  coefficients),  one  might  achieve  the 
desired  objective  by  fine  tuning  a  number  of  model  parameters  (e.g.,  the  order 


of  approximation  of  the  velocity  and /or  vorticity  gradients,  particularly  near 
the  wall,  mesh  size,  time  step,  type  of  discretization,  outer  boundary,  just  to 
name  a  few).  If  one's  objectives  are  to  perform  numerical  experiments  for 
sake  of  numerical  experiments,  with  no  concern  with  the  compatibility  of  the 
numerical  and  experimental  results,  then  one  can  objectively  asses  the  model 
instead  of  attempting  to  attribute  to  it  artificial  powers  of  prediction. 

As  far  as  the  turbulent  flows  are  concerned,  some  or  all  of  the  predictions 
of  the  numerical  calculations  for  a  given  flow  depend  on  the  closure  model 
used.  Some  models  do  better  than  others  for  some  flows  and  worse  than 
others  for  other  flows.  No  model,  however  sophisticated,  has  a  corner  on  the 
numerical  market.  Among  the  numerous  theoretical,  numerical,  and 
experimental  investigations,  impulsively-started  steady  flow  about  a  circular 
cylinder  has  occupied  a  prominent  place  partly  because  of  its  intrinsic  interest 
towards  the  understanding  of  the  evolution  of  separation,  vortex  formation, 
growth,  and  partly  because  it  provided  the  most  fundamental  case  for  the 
comparison  and  validation  of  various  numerical  methods  and  codes.  In 
recent  years,  attention  has  turned  to  a  broader  class  of  relatively  manageable 
time-dependent  flows  about  bluff  bodies:  Non-impulsively-started  flows, 
sinusoidally  oscillating  flows,  co-existing  flows  (uniform  flow  plus  oscillating 
flow),  flow  from  one  steady  state  to  another  (at  a  lower  or  higher  Reynolds 
number  through  the  use  of  prescribed  changes  in  velocity),  and  so  on.  The 
solution  of  these  problems  at  sufficiently  high  Reynolds  numbers  will  have 
far  reaching  theoretical  and  practical  consequences.  As  noted  above,  this  is 
not  yet  the  case,  and  the  solutions  must  necessarily  be  confined  to  cases  where 
the  accurate  prediction  of  physical  experiments  and  the  instant  gratification 


are  not  the  real  objectives.  However,  it  is  hoped  that  even  the  approximate 
solutions  will  have  enough  information  to  elucidate  the  physics  of  the 
phenomenon. 

The  hydrodynamic  loading  situations  which  are  well  understood  are 
those  which  do  not  involve  flow  separation.  Thus,  they  are  amenable  to 
nearly  exact  analytical  treatment.  These  concern  primarily  the  determination 
of  the  fluid  forces  on  large  objects  in  the  diffraction  regime  where  the 
characteristic  dimension  of  the  body  relative  to  the  wave  length  is  larger  than 
about  0.2.  The  use  of  various  numerical  techniques  is  sufficient  to  predict 
accurately  the  forces  and  moments  acting  on  the  body,  provided  that  the 
viscous  effects  and  the  effects  of  separation  for  bodies  with  sharp  edges  are 
ignored  as  secondary. 

The  understanding  of  the  fluid-structure  interactions  which  involve 
extensive  flow  separation  and  dependence  on  numerous  parameters  such  as 
Reynolds  number,  Keulegan-Carpenter  number,  relative  roughness,  relative 
motion  of  the  body,  proximity  effects,  hydroelastic  response,  etc.  is  far  from 
complete  (Sarpkaya  &  Isaacson  1981).  There  are  several  reasons  for  this.  First, 
although  the  physical  laws  governing  the  motion  (the  Navier-Stokes 
equations)  are  well  understood,  valid  approximations  necessary  for 
numerical  and  physical  model  studies  are  still  unknown.  Even  the 
unidirectional  steady  flow  about  a  bluff  body  remains  theoretically 
unresolved.  Much  of  our  understanding  of  vortex  shedding  behind  bluff 
bodies  came  from  steady-flow  experiments,  highly  idealized  models,  and 
limited  numerical  solutions.  Most  of  the  numerical  studies  based  on  the  use 
of  the  Navier-Stokes  equations  and  some  suitable  spatial  and  temporal 


differencing  schemes  are  limited,  out  of  necessity,  to  low  Reynolds  number 
flows.  A  second  reason  why  progress  has  been  slow  is  that  the  bluff  body 
problems  involving  wake  return  are  an  order  of  magnitude  more  complex 
and  there  has  been  only  a  handful  of  limited  applications  of  the  methods 
based  on  Navier-Stokes  equations. 

The  formation  of  a  wake  gives  rise  not  only  to  a  form  drag,  as  it  would  be 
the  case  if  the  motion  were  steady,  but  also  to  significant  changes  in  the 
inertial  forces.  The  velocity-dependent  form  drag  is  not  the  same  as  that  for 
the  steady  flow  of  a  viscous  fluid,  and  the  acceleration-dependent  inertial 
resistance  is  not  the  same  as  that  for  an  unseparated  unsteady  flow  of  an 
in  viscid  fluid.  In  other  words,  the  drag  and  inertial  forces  are  interdependent 
as  well  as  time-dependent.  These  effects  are  further  compounded  by  the 
diffusion  and  decay  of  vortices  and  by  the  three-dimensional  nature  of 
vorticity  due  to  turbulent  mixing,  finite  spanwise  coherence,  and  the  random 
nature  of  the  vortices  (which  give  rise  to  cycle-to-cycle  variations  and 
numerous  flow  modes  even  under  controlled  laboratory  conditions).  The 
stronger  and  better  correlated  the  returning  vortices,  the  sharper  and  more 
pronounced  the  changes  are  in  pressure  distribution  on  the  body  and  in  the 
integrated  quantities  such  as  the  lift,  drag,  and  inertia  coefficients. 

It  is  clear  from  the  foregoing  that  the  objectives  of  the  present 
investigation  are  to  carry  out  extensive  numerical  experiments  through  the 
use  of  the  vorticity-stream  function  form  of  the  Navier-Stokes  equations  and 
their  finite  difference  form,  on  co-existing  flows  (sinusoidal  oscillation  plus 
steady  mean  flow).  The  expectations  are  that  the  results  will  point  out  the 
strengths  and  weaknesses  of  the  code,  for  the  particular  type  of  formulation 


used,  explain  the  reasons  between  the  various  numerical  predictions  of  the 
same  problem,  and,  hopefully,  shed  some  light  on  the  physics  of  flows 
heretofore  uncalculated. 


II.    BACKGROUND  STUDIES 

A  finite  difference  analysis  of  the  Navier-Stokes  equations  for  a 
sinusoidally-oscillating  ambient  flow  about  a  circular  cylinder  at  K  (Keulegan- 
Carpenter  Number)  =  UmT/D  =  5  (Re  =  1000)  and  K  =  7  (Re  =  700)  has  been 
attempted  by  Baba  &  Miyata  (1987).  Their  results  have  shown  that  the 
calculations  can  be  carried  out  only  for  short  times  (less  than  two  cycles  of 
flow  oscillation)  with  a  non-super  computer.  Murashige,  Hinatsu  and 
Kinoshita  (1989)  have  used  a  similar  method  to  analyze  three  cases  (K  =  5,  7, 
and  10)  at  higher  Reynolds  numbers  around  104.  The  flow  was  perturbed  by 
artificial  means  to  trigger  an  asymmetry.  At  K  =  10,  a  transverse  vortex  street 
appeared,  in  agreement  with  experimental  observations.  The  numerical 
simulation  of  steady  flow  past  a  circular  cylinder  undergoing  in-line  and/or 
transverse  oscillations  through  the  use  of  two-dimensional  unsteady  Navier- 
Stokes  equations  was  undertaken  by  Lecointe  et  al.  (1987)  for  relatively  small 
amplitudes  (A/D  =  0.13).  Justesen  (1991)  presented  extensive  results  obtained 
from  a  numerical  solution  of  a  vorticity-stream  function  formulation  of  the 
Navier-Stokes  equations  for  the  flow  around  a  circular  cylinder  in  planar 
oscillating  flow  at  small  Keulegan-Carpenter  numbers  in  the  subcritical 
Reynolds  number  range.  Justesen  introduced  a  straining  parameter  "a"  in 
order  to  better  resolve  the  large  gradients  near  the  cylinder  surface.  This  is  in 
addition  to  the  logarithmic  straining,  commonly  used  as  part  of  the 
transformations,  for  a  better  resolution  of  the  gradients  near  the  body. 


Evidently,  Justesen's  transformation  for  a  =  0  defaults  to  the  logarithmic 
straining.  However,  "a"  becomes  another  disposable  parameter,  dependent 
on  at  least  K  and  Re.  Justesen  had  to  choose  judiciously  the  value  of  the 
straining  parameter  for  each  K  in  order  achieve  drag  and  inertia  coefficients 
in  satisfactory  agreement  with  those  obtained  experimentally.  A  systematic 
numerical  variation  of  the  governing  parameters  for  an  arbitrary  U(t)  is 
extremely  difficult. 

The  in-line  oscillations  of  a  cylinder  in  uniform  flow  (or  the  sinusoidally 
oscillating  flow  with  a  steady  mean  flow)  has  been  the  subject  of  intense 
interest  in  recent  years  (see,  e.g.,  Sarpkaya  &  Isaacson,  1981  and  Sarpkaya  & 
Storm,  1985)  in  connection  with  the  understanding  of  the  behavior  of  hot- 
wire anemometers  and  the  fluid  loading  of  structures  subjected  to  currents, 
gusts  and  other  types  of  unsteady  flows.  The  biassing  of  the  shedding  of  the 
vortices  by  the  current  causes  profound  changes  in  both  the  drag  and  inertia 
coefficients,  relative  to  their  no-current  values.  The  mobile  separation  points 
undergo  large  excursions,  as  much  as  120  degrees  during  a  given  cycle  of 
oscillation  over  a  circular  cylinder  (Sarpkaya  and  Butterworth,  1992).  These 
effects  are  further  compounded  by  the  diffusion  and  decay  of  vortices  and  by 
the  three-dimensional  nature  of  vorticity  due  to  turbulent  mixing,  reduced 
spanwise  coherence,  mutual-induction  instability,  and  the  random  nature  of 
vortices  which  give  rise  to  cycle-to-cycle  variations  and  numerous  flow 
modes  even  under  controlled  laboratory  conditions.  It  is  because  of  these 
reasons  that  the  present  work  is  dedicated  to  the  numerical  analysis  of  such 
flows. 


III.  NUMERICAL  REPRESENTATION 

A.    COMPUTATIONAL  METHOD 

Here  only  a  brief  description  of  the  computational  method  is  presented. 
A  more  in  depth  description  is  given  by  Wang  (1989)  and  Putzig  (1991). 

The  fluid  is  assumed  to  be  two-dimensional,  incompressible  and  viscous. 
The  governing  equations  for  the  solution  are  the  Navier-Stokes  equations, 
with  the  stream  function  and  the  vorticity  as  independent  variables.  To 
achieve  a  higher  density  of  mesh  points  near  the  cylinder  surface,  the 
computational  domain  is  transformed  from  the  physical  plane  in  Fig.  1  to  a 
rectangular  plane  in  Fig.  2  (see  Appendix).  In  the  rectangular  plane,  the  mesh 
is  maintained  at  a  uniform  grid  spacing.  It  is  necessary  to  have  more  mesh 
points  closer  to  the  cylinder  surface  because  in  this  region  the  gradients  of 
both  the  vorticity  and  the  stream  function  are  the  largest. 

A  third-order  in  time,  second-order  in  space,  three-level  predictor- 
corrector  finite-difference  scheme  is  used  to  solve  the  vorticity-transport 
equation.  A  Fast  Poisson  Solver  automatically  discretizes  the  separable 
elliptic  equation  which  is  then  solved  by  a  generalized  cyclic  reduction 
algorithm.  One  of  the  solvers  was  provided  by  the  FMSL  mathematics  library, 
for  the  solution  of  separable  elliptic  partial  differential  equations  with  a 
square  grid,  and  the  other  by  the  National  Center  for  Atmospheric  Research, 
for  the  solution  of  separable  elliptic  partial  differential  equations  with  a  non- 
square  grid. 
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The  unsteady  Navier-Stokes  equations  in  the  polar  coordinates,  as 
defined  by  the  vorticity  transport  equation  and  the  vorticity/ stream-function 
equation  are, 


3co_l 

at    r 


—  (CO— -) (CD— -) 

dr     ae     ae     ar 


=  vV2co 


(1) 


and 


where 


V2v  =  co  (2) 


n2    a2    i  a     i  a2 
v  =  — =•+  — +■ 


ar2    rar   r2ae2  (3) 

co  and  \j/  are  the  vorticity  and  the  stream  function,  v  is  the  kinematic  viscosity, 
t  is  the  time  and,  r  and  0  are  polar  coordinates  directions  (see  Fig.  1).  The 
velocity  components  in  the  r  and  0  directions  are  defined  by 

1 3\i/  3w 

u  = —  v  =  — - 

r  30     and  3r  (4) 

The  boundary  conditions  for  the  physical  problem  are: 
(1)  no  slip  and  zero  normal  velocity  on  the  surface  of  the  cylinder 

\i/  =  -^  =  0     on     r  =  R 

*  (5) 

and  (2)  the  potential  flow  at  infinity  is  defined  as 


w  =  U(r )sin0 

r  (6) 

and  co  =  0  at  r  =  co.  u  is  the  external  flow  and  R  is  the  radius  of  the  cylinder. 


The  coordinate  transformations  required  to  go  from  the  physical  domain 
to  the  computational  domain  are: 

r  =  R  exp(a^)     and     6  =  arj  (7) 

where  R  is  the  radius  of  the  cylinder  and  'a'  is  a  transformation  parameter. 

The  transformation  of  the  non-dimensionalized  vorticity-stream 
function  equations  and  their  finite  difference  form  through  the  use  of  the 
central  difference  approximation  for  vorticity  and  a  two-step,  three-level, 
predictor-corrector  scheme,  with  a  third  order  accuracy  in  time,  are  described 
in  detail  in  Wang  (1989),  in  Fredrickson  (1990)  and  in  Putzig  (1991)  and  will 
not  be  repeated  here. 

B.    CALCULATION  OF  THE  FORCE  COEFFICIENTS 

The  in-line  and  transverse  force  coefficients  are  determined  from  the 
combined  contributions  of  the  shear  and  pressure  forces  acting  on  the 
cylinder.  The  viscous  forces  are  calculated  from  ts=  fico.  The  total  in-line 
force  then  reduces  to 


fil  =  -JoKPs  cos(8)Rd6  -  /oKM.cosin(co)Rde  (g) 


and  the  total  lift  force  as 

Fl  =  -J02npssin(e)Rde-  j027tuxocos(co)Rde  (9) 

After  dividing  the  in-line  and  the  lift-force  equations  by  (0.5  pU2  D)  and  defining 
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-   _2(ps-Poo) 
Ps  "  tt2 


(10) 


the  force  coefficients  reduce  to 


CIL  =  -ij027Cpscos(6)Rde-—  J02n^a5sin(co)Rde 
2  Re 


(11) 


and 


CL=-ij02npssin(6)Rde-  — j02npLC0cos((0)Rde 


(12) 


The  pressure  coefficient  is  determined  from  the  Navier-Stokes 
equations  in  terms  of  dimensionless  vorticity.  Once  integrated  with  respect  to 
0,  one  has 


ps<e)  =  ps(0)+A£?*(^l     de 


ReJ0    VarJr=1 


(13) 


Equation  (13)  is  substituted  into  equations  (11)  and  (12)  to  determine  the 
numerical  scheme  for  the  in-line  and  transverse  force  coefficients, 


CIL  = 


2       f2K 


Re  Jo 


\\\ 

Jo        J( 


Ydco^ 


H  3r 


de 


A=i 


cos(8)  +  G)sin(8)ki8 


(14) 


and 


Re 


f{ 


de 


A=i 


sin(e)  +  wcos(e)^de 


(15) 
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The  radial  derivative  of  the  vorticity  on  the  surface  of  the  cylinder, 
appearing  in  Equations  (14)  and  (15),  is  determined  through  the  use  of  discrete 
pointwise  approximations  of  various  orders,  ranging  from  second  to  tenth 
order.   For  a  second  order  approximation,  one  has 


^tQj"]       =  -3o)i  +  4coi+1  -coi+2  ,  0(^2) 
^  3r  Jr=1  2Ac, 


(16) 


For  higher  order  polynomials  Equation  (16)  may  be  written  as 


fay 


Aco,  +  Bco1+1  +  Ccoi+2  +  Dcoi+3  +  Eco1+4  +  . 


At 


+  0(A^n)     (17) 


in  which  the  coefficients  A  -K  are  given  in  Table  1  below. 
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Table  1:  Coefficients  of  the  Polynomial  in  Eq.  (17) 

n=2  n=4  n=6  n=8  n=10 


A 

-3/2 

-25/12 

-49/20 

-761/280 

-7381/2520 

B 

2 

4 

6 

8 

10 

C 

-1/2 

-3 

-15/2 

-14 

-45/2 

D 

4/3 

20/3 

56/3 

40 

E 

-1/4 

-15/4 

-35/2 

-105/2 

F 

6/5 

56/5 

25^5 

G 

-1/6 

-14/3 

-35 

H 

8/7 

120/7 

I 

-1/8 

-45/8 

J 

10/9 

K 

-1/10 
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C     CALCULATION  OF  THE  DRAG  AND  INERTIA  COEFFICIENTS 

If  one  were  to  associate  the  total  force  with  a  velocity-square-dependent 
drag  force  and  an  acceleration-dependent  inertial  force  then  the  coefficient 
associated  with  the  latter  may  be  interpreted  as  some  measure  of  the  added 
mass.  But  one  must  bare  in  mind  that  such  a  decomposition  is  far  from  being 
unique. 

It  has  been  customary  to  express  the  fluid  force  acting  on  a  body  moving 
in  a  fluid  otherwise  at  rest  as 

F(t)  =  ipCS  Ap|{(U0  +  U(t)}|{(U0  +  U(t)}  +pkH  V^jp  (18) 

where  Uo  represents  the  steady  velocity;  U(t),  the  time-dependent 
oscillations;    Cj  ,  the  Fourier-averaged  drag  coefficient  and  k^  ,  the  Fourier- 

averaged  added-mass  coefficient.    It  is  customary  to  use  an  inertia  coefficient 
Cm  for  a  fluid  in  motion  about  a  body  at  rest  through  the  use  of  Cm  =  1  +  kV. . 

The  Fourier  averages  of  the  drag  and  added-mass  coefficients  over  a 
period  of  T  may  be  calculated  by  multiplying  both  sides  of  Equation  (8)  once 
with  U(t)  and  once  with  dU/dt  to  yield 

T 
2jF(t)U(t)dt 

Q  =  — T s (19) 

pApJ|{(U0  +U(t)}|{(U0  +U(t)}U(t)dt 
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and 

T 


dUft) 
dt 


fF(t)^dt 


(Cm  -  1)  =  k^  =  — ^ (20) 


pvJ[dU(t)/dt]dt 
0 


which  may  be  evaluated  readily  provided  that  sufficiently  reliable  data  are 
available  for  F(t),  U0,  U(t),  and  dU(t)/dt. 

A  simple  dimensional  analysis  of  the  flow  under  consideration 
shows  that  the  time-averaged  force  coefficients  (CjJ  and  kV1-  )  are  functions  of 

a  relative  amplitude  or  Keulegan-Carpenter  number,  Mach  number, 
Reynolds  number,  and  a  parameter  involving  Uo  (e.g.,  UoT/D  or 
Uo/[U(t)]max)-  There  are  numerous  possibilities  regarding  the  definitions  of 
the  relative  amplitude  or  Keulegan-Carpenter  number  and  the  Reynolds 
number.  The  purpose  of  the  search  for  a  more  suitable  Keulegan-Carpenter 
number  and /or  Reynolds  number  is  to  enhance  the  correlation  of  the  data  to 
reduce  the  number  of  the  governing  parameters,  possibly  eliminating  UoT/D 
as  an  independent  parameter.  The  list  of  possible  Reynolds  numbers  and 
Keulegan-Carpenter  numbers  is  long  and  will  not  be  given  here.  Suffice  it  to 
note  that  the  two  force-coefficients  for  the  flow  about  a  cylinder  may  be 
written  as 

Cu 

ud    =fi(K/Re/VK)  (21) 

kij 
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or  as 


=  fi(K+,Re+,VK)  (22) 


in  which 

K  =  UmT/D     ,     Re  =  UmD/v     ,     VK  =  U0T/D 
K+  =  K(l  +  Uo/Um)  =  K  +  U0T/D 
Re+  =  Red  +  Uo/Um)  =  Re  +  U0D/v  (23) 

The  purpose  of  the  present  calculations  was  not  to  provide  a  detailed 
comparison  between  the  measured  and  calculated  forces  but  rather  to  attempt 
to  establish  a  relationship  between  the  shedding  of  vortices  and  the  relative 
magnitude  of  the  current.  The  particular  values  of  K  and  Vr  chosen  for  the 
calculations  (K  =  4  -  6,  Vr  =  0.0  -  1.2)  was  one  for  which  some  experimental 
data  were  available  at  comparable  p,  Re,  and  Vr  values. 
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IV.   DISCUSSION  OF  RESULTS 

The  numerical  experiments  were  carried  out  through  the  use  of  a  VAX- 
2000,  a  VAX-3520,  a  CRAY  Supercomputer.  The  solution  procedure  and 
technique  have  been  validated  for  several  types  of  unsteady  flows,  i.e., 
impulsively-started,  suddenly-stopped,  and  uniformly-decelerated  flows 
before  applying  it  to  co-existing  flows  (oscillatory  flow  or  sinusoidally- 
oscillating  flow  with  mean  velocity).  Excellent  agreement  with  flow 
visualization  and  experimentally  determined  drag  and  lift  coefficients  has 
been  obtained  for  both  symmetric  and  asymmetric  wake  solutions.  It  is  this 
validation  that  led  to  the  exploration  of  the  characteristics  of  sinusoidally- 
oscillating  flows  superimposed  on  a  mean  velocity.  The  oscillation  was 
specified  by  U  =  U0+Um  sin(27ct/T)  in  which  U0  is  the  steady  mean  velocity 
and  Um  is  the  amplitude  of  sinusoidal  oscillations. 

The  flow  was  perturbed  by  changing  the  direction  of  the  ambient  flow 
sinusoidally  (with  an  amplitude  of  one-half  of  a  degree)  during  the  first  cycle 
of  the  oscillation.  The  amplitude  of  the  sine  wave  was  the  only  free 
parameter.  It  is  worth  noting  that  this  type  of  disturbance  gradually  returns 
the  perturbed  quantity  to  its  initial  state. 

Numerical  experiments  have  been  carried  out  in  the  range  of  K  =  1-4,  (3 
=  200,  Re  =  200-800,  At  =  0.00025-0.0005  and  for  various  values  of  Vr  =  U0/Um, 
as  shown  in  Table  2. 
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Table  2  The  range  of  the  governing  parameters 

(P  =  200) 


Vr    / 

K=  1 

2 

2 

4 

0.0 

X 

X 

X 

X 

0.60 

X 

0.65 

X 

0.70 

X 

Evidently,  this  is  a  rather  limited  exploration  of  a  highly  complex 
problem  and  requires  much  more  numerical  and  experimental  work.  The 
purpose  of  the  present  calculations  was  not  to  provide  a  detailed  comparison 
between  the  measured  and  calculated  forces  but  rather  to  attempt  to  establish 
a  relationship  between  the  shedding  of  vortices  and  the  relative  magnitude  of 
the  current. 

Figures  3  through  6  show  for  K  =  1,  £  =  200,  Vr  =  0.0,  and  At  =  0.00025,  the 
in-line  force  and  streamlines.  The  calculated  force  coefficients  are  Cd  =  2.18 
and  Cm  =  2.15.  The  corresponding  theoretical  values  are  Cd  =  2.2  and  Cm  =  2.1 
(Wang  1968).  Figures  7-13  show  the  in-line  force,  the  streamlines  and  the 
streaklines,  for  K  =  2,  p  =  200,  and  Vr  =  0.0.  The  calculated  force  coefficients 
are  Cd  =  1-35  and  Cm  =  2.13.  The  corresponding  theoretical  and  experimental 
values  are,  respectively,  Cd  (th)  =  0.95,  Cd  (exp)  =  1.50  and  Cm  (th)  =  2.02  and  Cm 
(exp)  =  2.0.  The  differences  between  the  calculated,  experimental,  and 
theoretical  values  are  due  to  the  fact  that  the  theoretical  values  do  not 
account  for  the  flow  separation.    Figures  14-24  show  the  in-line  force  and 
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detailed  streakline  plots  for  K  =  3,  (3  =  200,  and  Vr  =  0.0.  The  calculated  force 
coefficients  are  Cd  =  1.35  and  Cm  =  2.02.  The  corresponding  experimental 
values  (Bearman  et  al.  1985)  are  Cd  =  1.45  and  Cm  =  2.10.  There  are  no 
theoretical  values  to  compare  with  for  the  reasons  just  cited.  Figures  25-31 
show  the  in-line  force  and  the  streaklines  for  K  =  4,  p  =  200,  and  Vr  =  0.0.  The 
calculated  force  coefficients  are  Cd  =  1-38  and  Cm  =  1.92.  The  corresponding 
experimental  values  (Bearman  et  al.  1985)  are  Cd  =  1-45  and  Cm  =  1.95.  It  is 
clear  from  the  foregoing  that  for  the  no-current  case  the  calculated  force- 
transfer  coefficients  are  in  very  good  agreement  with  those  obtained 
experimentally.  For  the  cases  for  which  a  comparison  is  possible,  the 
calculated  values  also  agree  reasonably  well  with  the  theoretical  values  of 
Wang  (1968). 

Figures  32-38  show  the  in-line  and  transverse  forces,  the  streamlines  and 
the  streaklines  for  K  =  4  and  Vr  =  0.6  (with  Cd=  1.17,  Cm  =  1.65);  and  Figs.  39-41 
show  the  in-line  and  transverse  forces  and  a  single  streakline  plot  for  K  =  4 
and  Vr  =  0.65  (with  Cd=  1.16,  Cm  =  1.63);  and  finally,  Figs.  42-43  show  a 
streamline  and  a  streakline  plot  for  K  =  4  and  Vr  =  0.7  (with  Cd=  1.15,  Cm  = 
1.63). 

The  calculation  with  current  were  carried  out  for  the  expressed  purpose  of 
substantiating  the  previous  findings  that  for  K  =  4,  within  a  narrow  range  of 
Vr  values,  the  width  of  the  wake  increases  and  the  vortices  begin  to  arrange 
themselves  along  three  rows.  In  order  to  make  sure  that  this  finding  was  not 
a  consequence  of  the  sensitivity  of  the  code  to  grid  shape  and  size  and  the 
time-interval   used,   a   more   comprehensive   series   of  calculations   were 
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undertaken  through  the  use  of  a  non-square  grid.  Not  only  were  the 
previous  conclusions  substantiated  but  also  additional  facts  were  uncovered. 
A  comparison  of  Fig.  38  (for  Vr  =  0.60),  with  Fig.  41  (for  Vr  =  0.65),  and  Fig.  43 
(for  Vr  =0.7)  shows  that  the  wake  is  comprised  of  three  rows  of  heterostrophic 
vortices.  They  differ  only  in  detail  from  one  Vr  to  another  in  the  narrow 
range  of  Vr  values  from  about  0.6  to  0.7.  At  lower  Vr  values,  the  inner  pair  of 
vortices  are  closer  to  the  wake  axis.  As  Vr  increases  toward  0.7,  the  inner  pair 
reaches  the  edges  of  the  outer  pair  (Fig.  43),  and  the  inner  and  the  outer 
vortices  form  an  interesting  quadruple.  The  vortex  pairs  on  one  side  of  the 
street  may  propel  themselves  (as  a  Kelvin  oval)  away  from  the  wake  axis,  or 
become  part  of  a  larger  scale  instability.  In  either  case,  it  will  be  difficult  to 
predict  the  behavior  of  such  a  special  wake  at  larger  times  because  it  is  rather 
difficult  to  distinguish  between  the  naturally  occurring  instabilities  and  those 
occurring  numerically  due  to  truncation  errors.  Nevertheless,  it  is  gratifying 
to  note  that  the  flow  visualization  of  the  early  stages  of  the  wake  for  K  =  4,  at  a 
representative  value  of  Vr  =  0.65,  yielded  results  in  excellent  agreement  with 
those  presented  herein. 

Table  3  shows  a  comparison  of  the  calculated  and  experimental  (Moe  & 
Verley  1980)  drag  and  inertia  coefficients  for  K  =  4  for  representative  values  of 
Vr.  As  expected,  the  inertia  coefficients  agree  extremely  well.  As  far  as  the 
drag  coefficients  are  concerned,  the  agreement  is  not  as  good  but  certainly 
better  than  expected  in  view  of  the  fact  that  (3  was  200  in  the  calculations  and 
about  600  in  the  experiments.  Nevertheless,  the  trend  of  the  data  is  well 
predicted. 
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Table  3.        Comparison  of  the  Measured  and  Calculated  Force 
Coefficients  for  K  =  4  and  Three  Values  of  Vr 


Vr: 

0.60 

0.65 

0.70 

Cm  (exP> 

1.68 

1.63 

1.60 

C"   (cal) 
m 

1.65 

1.63 

1.63 

Cjj  (exp) 

0.93 

0.98 

1.05 

C"  (cal) 

1.17 

1.16 

1.15 
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V.  CONCLUSIONS 

The  investigation  reported  here  warranted  the  following  conclusions: 

1.  Even  the  higher  order  finite  difference  formulations  of  the  governing 
equations  based  on  the  vorticity/stream-function  formulation  of  the 
Navier  Stokes  equations  can  be  solved  for  only  relatively  small 
Reynolds  numbers.  This  is  primarily  due  to  stability  and  computer 
constraints. 


The  numerical  experiments  with  pulsating  flows  (oscillation  plus 
steady  mean  flow)  for  K  =  4  yielded  force-transfer  coefficients  in  good 
agreement  with  those  obtained  experimentally. 


For  K  =  4  and  relative  current  velocities  of  about  0.6-0.7,  the  vortices 
shed  nearly  symmetrically  at  each  cycle  and  gave  rise  to  a  most 
unusual  three-row  vortex  street,  where  each  row  is  comprised  of  a 
pair  of  heterostrophic  vortices.  For  relative  current  velocities  larger 
than  about  one,  the  vortex  wake  returned  to  the  asymmetric  mode,  as 
is  encountered  in  a  regular  Karman  vortex  street.  The  use  of  a  non- 
square  grid  in  lieu  of  a  square  grid  in  solving  the  Poisson's  equation 
did  not  alter  the  results. 


The  foregoing  numerical  experiments  could  not  have  been  possible 
had  it  not  been  due  to  the  availability  of  a  VAX-3520  and  a  CRAY 
supercomputer.  It  is  fully  realized  that  numerical  instabilities  versus 
fluid  dynamical  instabilities  have  different  and  at  times  competing 
mathematical  and  computational  demands.  It  is  because  of  this 
reason  that  calculations  at  higher  Reynolds  numbers  (while 
preserving  flow  stability)  and  calculations  at  very  small  Keulegan- 
Carpenter  numbers  (while  preserving  numerical  stability)  require 
extremely  large  CPU  times  even  on  a  supercomputer. 
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Figure  1.    Grid  in  the  Computational  Domain 
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Figure  2.  Grid  in  the  Physical  Domain 
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Figure  3.  In-Line  Force  Coefficient,  K  =  1,  Re  =  200,  Vr  =  0.0 
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Figure  7.  In-Line  Force  Coefficient,  K  =  2,  Re  =400,  Vr  =  0.0 
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Figure  8.  Streamlines,  K  =  2,  Re  =400,  Vr  =  0.0,  t/T  =  4.0 
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Figure  9.  Streamlines,  K  =  2,  Re  =400,  Vr  =  0.0,  t/T  =  5.0 
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Figure  10.    Streamlines,  K  =  2,  Re  =400,  Vr  =  0.0,  t/T  =  6.0 
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Figure  14.    In-Line  Force  Coefficient,  K  =  3,  Re  =600,  Vr  =  0.0 


36 


o 

CO 


-»-> 

o 

o 

II 

t-, 

> 

o 
o 

II 

o; 

CO 

II 

to 
a; 

C 

OS 

a; 

to 

en 
in 

H 

O) 

u 
3 

be 

•  ^^ 

tin 


37 


o 

II 
H 

■*-> 

o 

o 

II 

> 

o 

o 

II 

to 
II 

en 

0> 

G 

•  ■-"« 

«J 
<U 

is 

CD 

vO 

T— I 

u 

3 
bO 

K 


38 


II 
H 


o 

© 

II 


o 

o 

II 

PS 

en 
II 

<u 

c 

<u 

CD 


bo 


39 


o 
in 

II 
H 


o 
o 


> 

o 
o 

II 

& 

CO 


in 


2 

0) 

is 

00 
t— < 

(V 
U 

3 

to 


40 


•  114 


IT) 
IN 

II 

H 


o 

o 

II 


o 
o 

II 

CO 


C 
CD 


t— i 
M 

to 

tin 


41 


o 

in 

ii 

H 

-*— > 

o 

o 

ll 


o 
o 

II 

& 

CO 

II 

<U 
C 

cu 

Cf) 


o 


42 


in 

II 

H 


o 
o 


> 

o 
o 


g 

■s 

en 


3 

fad 


43 


o 
in 


o 
o 


> 

o 

o 

II 

CO 


en 

c 
en 


(N 
.■-I 


44 


in 

ii 
H 


o 

o 

II 


o 

o 

II 
CO 

II 

cu 

C 

if) 


5Q 

be 


45 


o 

II 
H 

o 

o 

II 


o 

o 

II 

o 

CO 

II 

I/} 

a 
C 

a 
is 

CO 


Hi 

to 

•  *■* 


46 


\ 

CO 

u!  o 

i  A 

O  • 

O  * 

UJ 

\ 

INUNE  F< 

763        -0.5 

1 

N 

• 

01 

o 

• 

V 

0.0005 


1.2504 


2.6003 

t/T 


3.7501 


Figure  25.    Inline  Force  Coefficient,  K  =  4,  Re  =  800,  Vr  =  0.0 
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Figure  32.    Inline  Force  Coefficient,  K  =  4,  Re  =  800,  Vr  =  0.6 
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Figure  33.    Transverse  Force  Coefficient,  K  =  4,  Re  =  800,  Vr  =  0.6 
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Figure  39.    Inline  Force  Coefficient,  K  =  4,  Re  =  800,  Vr  =  0.65 
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Figure  40.    Transverse  Force  Coefficient,  K  =  4,  Re  =  800,  Vr  =  0.65 
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